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Extensive simulations are performed of the diffusion-limited reaction A+B— > in one dimension, 
with initially separated reagents. The reaction rate profile, and the probability distributions of 
the separation and midpoint of the nearest-neighbour pair of A and B particles, are all shown to 
exhibit dynamic scaling, independently of the presence of fluctuations in the initial state and of an 
exclusion principle in the model. The data is consistent with all lengthscales behaving as f^''* as 
t ^ 00. Evidence of multiscaling, found by other authors, is discussed in the light of these findings. 

PACS numbers: 05.40. -fj, 02.50.-r, 82.20.-w 



I. INTRODUCTION 

There has been a lot of recent interest in the scal- 
ing behaviour of the reaction front that exists be- 
tween regions of initially separated reagents A and B 
that perform Brownian motion and annihilate upon 
contact according to the reaction scheme A-|-B^ 
§|li,||,0|.|.0,0E|jll,0,g,^,p. The evolu- 
tion of the particle densities ajx^t) and bjx, t) at position 
X and time t is governed by the equations 



da d^a 



db _ j^dH 



at 



R, 



where R is the reaction rate per unit volume, and the 
diffusion constant D has been assumed equal for both 
species. For the Boltzmann equation ansatz R — kab, 
the solution to the resulting partial differential equations 
with the initial condition 



a{~x,0) — ao9{x) = b{x,0), 



(2) 



where 9 is the Heavyside function, has the scaling prop- 
erty 



^^^-a-(l/2)^|^|.^ forX«il/2, 



(3) 



with a = i Q . This resuh may be understood by consid- 
ering the steady-state solutions to Eqs. ([^) for boundary 
conditions a{—x) — > J\x\/D, a(x) —^ 0, b{—x) — > and 
b{x) — > Jx/D 'AS X 00. Under these conditions, there 
are opposing constant currents J of either species, and it 
can be shown ||^,|l^ that the resulting reaction profile is 
of the form 



with A = 3 Returning to the the time-dependent case, 
the quantity (a — b) obeys a diffusion equation, whose 
solution for initial conditions (0) is 
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Let us assume that the reaction takes place within a re- 
gion of width w ^ t", with a < ^. The profiles for 
w X t^/"^ arc of the form a oc aox/t^^^, so there is 
a current of particles arriving the origin of the form J = 
Ddxtt ~ t~^^^. The characteristic timescale on which this 
current varies is (dlog J/dt)^^ cx t, whereas the equilibra- 
(1) tion time for the front is of order (w'^/D) '-^ t^" <C t. The 
front is therefore formed quasistatically, and so Eq. (H) 
may be obtained from ^ simply by writing J oc . 

Simulations and experiments appear to confirm these 
results when the spatial dimension d is two or greater 
lll'ilJlJll- dimension less than two, strong correlations 
between the motions of the two species cause the Boltz- 
mann approximation R ~ kab to break down. However, 
the solution to the steady state problem is still of the 
~|), albeit with a different exponent A — l/{d+ 1) 
If the results from the steady-state may still be 
used, this would lead again to dynamical scaling of the 
form (^, with a = j in d — 1. Simulations using a 
one-dimensional Probabilistic Cellular Automata (PCA) 
model appeared to verify the dynamical scaling form (|^) , 
though with a = 0.293 ± 0.005 (|L Monte-Carlo simula- 
tions also found a « 0.30 ± 0.01 [§. 

However, a recent article by Araujo et al ||l^] has chal- 
lenged the validity of the scaling form (||) . This article 
reported Monte Carlo (MC) simulations in one dimen- 
sion, using an algorithm where the A and B particles 
always react on contact and so are unable to cross over 
each other. The right-most A particle (RMA) is there- 
fore always to the left of the left-most B particle (LMB). 
Defining Iab as the separation between the RMA and 
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LMB, and m as the midpoint between them, Araujo et 
al found that the probabihty distributions Pi and Pm 
of these quantities displayed dynamic scaling, with char- 
acteristic lengthscales i^^^ and t^^^ respectively. Mean- 
while, the different moments of the reaction profile were 
described by a continuous spectrum of lengthscales be- 
tween i^/^ and t^^^. More specifically, defining 
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R{x, t)dx 
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Araujo et al found that Z^''-' 



i"" with 4 < a„ < 



- ^(q) _ ^3/8^ but 

I increasing monotonicaUy 



with q. They also proposed the following form for R: 



R{x,t)^t~'^/* 



exp 
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which predicted values of aq that were in good agreement 
with their numerical findings (the prefactor essen- 
tial for consistency, is missing in [p^). The authors of 
I p^ argued that Poisson noise in the initial state causes 
the reaction centre to wander anomalously as to cx 
invalidating the use of the steady state results. 

In this paper, I first describe extensive simulations of 
this system, using two independent models — the PCA 
model used in H, and a MC model similar to that of 
Araujo et al in ]16|. I find that dynamical scaling ap- 
pears to hold for Pi, Pm, and R, independently of the 
existence of Poisson fluctuations in the initial state and 
of the presence of an exclusion principle. While I con- 
firm the result Z''-' ^ t^/^, I find instead that both m'^''^ 
and x^*) appear to scale as t", with a w 0.28 ± 0.01, 
independently of q. The high statistics and wide time 
domain accessible in the PCA simulations show that this 
exponent is decreasing monotonicaUy in time, consistent 
with the asymptotic result a — ^ predicted by the anal- 
ogy with the steady-state result. The measured forms of 
Pi, Pm and R are found to be described by very simple 
analytic forms to high accuracy. I then discuss the va- 
lidity of the fluctuation argument used by Araujo et al 
to explain the result to ~ t^/^. An exact calculation of a 
related quantity suggests that the wandering of the reac- 
tion centre should instead be of the order ^ i^^^, which is 
not sufficient to make the use of the steady-state analogy 
invalid. Some of these results have been discussed in a 
previous publication jTsl. 



II. MONTE-CARLO MODEL 
A. Description of Model 

The model described in Ref. |l^ consists of indepen- 
dent random walkers with no exclusion principle. In 
the interests of computational efficiency, I used a model 
which is identical provided the site occupation number 
is not too large, but whose site updates may be effected 
using a lookup-table algorithm. In this way, it was pos- 
sible to obtain statistics equivalent to the simulations in 
]l6t in the space of a few days. 

The model has an 'exclusion principle', in that no more 
than 2lp particles of each type are allowed per site. In 
the diffusion step, each of these particles moves onto a 
neighbouring site, in such a way that no more than Ip 
particles may move from a given site in the same direc- 
tion at once. This constraint automatically satisfies the 
'exclusion principle'. If there are Ip or less particles on 
a site, then the direction in which each particle moves is 
chosen independently at random. If there are more than 
Ip particles, the same redistribution method is used for 
the 'holes' — i.e. the probability of j particles moving to 
the right when the occupation number is k is the same as 
[Ip — j) particles moving to the right when the occupa- 
tion number is {2lp ~ k). The diffusion constant for this 
model is ^. 

In these simulations, the value Ip — 13 was used (this 
was the largest value that could be implemented effi- 
ciently). Since the average density was 1 or less, the 
frequency of events where the occupation is greater than 
Ip is of order e~^/ {lp + l)\ « 4 x 10~^^, so these events are 
extremely rare (the simulations represent 21000 samples 
of 4000 sites over 25000 timesteps, so the expected total 
number of such events is less than 10). The influence 
of such events on the results is still smaller, since the 
probability of a large number of particles spontaneously 
moving in the same direction is low (e.g. 14 indepen- 
dent walkers move in the same direction with probabil- 
ity 2^^^^ Ri lO^''). Moreover, the universality class for 
the scaling properties is not expected to depend on such 
events, as the reactions take place in the zone where the 
density is low. These results may therefore justifiably 
described as equivalent to those reported in [Q. The 
Fortran implementation of this algorithm performed 
1.4 X 10'^ site updates per second on a HP 9000/715/75 
workstation. 

One timestep consists of moving all the particles, fol- 
lowed by a reaction step. The pure diffusion algorithm 
has a spurious invariance, in that particles initially on 
even sites will always be on even sites after an even num- 
ber of timesteps, and will be on odd sites at odd timesteps 
(and contrarily for particles initially on odd sites). In ac- 
cordance with the prescription in pq ] that an A particle 
never be found to the right of a B particle, it is important 
that the reaction takes account of particles of different 
types crossing over each other (i.e. an A-particle at site 
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FIG. 1. Illustration of the Monte Carlo reaction-diffusion 
algorithm, showing the diffusion step and the two stages of FIG. 2. Scaling plot of the density difference 0(2:, t) — 

the reaction step that first remove particles that have crossed for the MC data with random initial condition. Inset: bias 

over and then react those that are at the same site. plot for the total number of C particles Ec = / C{x, t) dx. 



i hopping to i + 1 at the same time that a B-particle 
at I + 1 hops to site i). The reaction algorithm first re- 
moved such particles, then annihilated any remaining A 
and B particles occupying the same site. This is illus- 
trated in Fig. |l|. Four sites are shown, with initially two 
sites occupied by A-particles (represented above the line, 
labeled 1-5) and two occupied by B-particles (beneath 
the line, labeled 6-10). In the diffusion step, each par- 
ticles moves onto a neighbour at random, producing the 
state (ii). The reaction first deletes the A-B pair (4,6) 
that crossed over, leading to state (iii), then removes the 
pairs that sit at the same site, giving a final state (iv). If 
the 'delete crossing' step were not present, the reaction 
step would simply remove the pair (4,8) from state (ii), 
leading to a state where there are B-particlcs to the left 
of A particles. 



B. Simulation Results 

An approximation to a Poissonian initial state of av- 
erage density unity was prepared by performing 16 at- 
tempts to add an A-particle, with probability 1/16, to 
each of the first 2000 sites of a 4000-site lattice. The 
other half of the lattice was similarly populated with B- 
particles. At the boundaries, particles that attempted to 
leave the system were allowed to do so, but a random 
number (distributed binomially between and 16, aver- 
age i) of particles was allowed to re-enter the system at 
the end sites. The average density at the extremities was 
thus kept at the value unity. 

In order to mimic the simulations in as closely 
as possible, instantaneous measurements were made of 
Iab^ ™, and the concentration profiles of the product 
and reagents at times 1000, 2500, 5000, 7500,. . . , 25000. 
These were then averaged over 21000 independent initial 
conditions. The quantities Z^'?-' and m^*^ were measured 
from the probability distributions over the samples, and 
a quantity A^^^ was defined as 
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where C = J Rdt is the profile of the reaction product. 
This quantity differs from x{q), but since J C dx oc t-^^^, 
and provided that x^"^-* behaves as a power of t, Eq. (5) 
of ||l^ shows that they should have the same scaling be- 
haviour. 

From Eq. (|^), the difference in the particle densities, 
(a — 6), obeys a simple diffusion equation, whose solution 
is given by (|^). Any finite-size effects in the data would 
first show up in deviations of the particle profiles from 
the values they would have for an infinite system. Figure 
Ijshows a plot of (a — b) as a function of (x/t^^^) for three 
time values, displaying excellent rescaling. A simpler test 
of finite-size effects is to show that the total C-particle 
number Sc = J C dx is proportional to t^^^. The inset 
to Fig. H confirms that Sc • t~^/'^ is indeed independent 
of time. 

Figure ^ is a log-log plot of X 
function of time, where 



and li^^ as a 



xi''\t) ^ C.x^'^Ht) 



(t)EE^X')(t) 



(11) 
(12) 

(13) 



and S.q, iJ-q, and Xq are constants that will be defined 
later. The straight lines are fits to the last 8 points for 
xi'^\ rrS^\ and The gradients for least squares fits 
to the curves in Fig. ^ are listed in Tabic |. The exponent 
describing Z^"^^ is close to as was found in ||l^. However, 
the results for tvP') and A^^^ differ dramatically from 
those of Araujo et al. Firstly, the exponent describing 
m{t) appears to be close to 0.29, instead of 0.375 as they 
found. Secondly, the exponents describing X^'^'^ appear 
to be independent of q. This means that C{x^t), and 
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FIG. 3. Log-log plot of ml"' 



and xi''^ (see text) from 
the MC simulations with Poisson initial conditions. 
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FIG. 4. Effective exponents for m*^', and X*^' (see 
text) from the MC simulations. 



by implication R{x,t), obeys a simple scaling form, in 
contrast to the anomalous form (||). 

To investigate for a trend in the exponents describing 
these quantities, the effective exponent (defined as the 
gradient between successive points in Fig. ||) is plotted 
as a function of l/log]^Q(f) in Fig. ^. The data for Z^^^ and 
rn'^' are far too noisy for any information to be obtained. 
The exponent for X'^^^ appears to to decrease slowly in 
time, but the time window in these simulations is too 
narrow for conclusive deductions to be made. 

Figures ||, ||, and are plots of C, Pmim) and Pi{Iab), 
as a function of appropriate scaling variables, to show 
the subjective quality of scaling for these quantities. 
The profiles of Pm{'m) and Pi{Iab) suggest the follow- 
ing forms: 
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FIG. 5. Scaling plot for C{x,t), for the MC simulations. 
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These forms predict the following results for the moments 
of these distributions: 



where 



1/9 



_ r /i9+i(V^/2)i/'? for g odd. 



(16) 

(17) 
(18) 



for q even. 

Using these values of fiq and A, in Eqs. (|ll|-[r^), one 
would expect ml''' and li'''^ to be independent of q if the 
forms ( [r^ , p^ ) are valid. The coincidence of the curves in 
Fig. 1^ confirms this. 

Figure ^ is an explicit test of the forms ( jliHTsI ) 
against the data, by plotting log[m(^)Pm(TO)] and 
log[Z(2)rip,(/)] against (m/m^^))^ and (l/l'-^yf respec- 
tively, at t = 25000. The Y-ordinate has been shifted so 
that all curves are coincident at the origin. The curve 
labeled Kmc is fog[C(x, 25000) - C(x, 22500)], which is 
approximately proportional to R{x, 25000), as a function 
of {x/X'-^^)'^. The straight line for this curve suggests 
that the reaction profile R{x,t) is also a Gaussian. This 
again contradicts the form (0) proposed by Araujo et al. 
It is not possible, however, to derive analytical forms for 
S^q that lead to xi''^ being independent of q without as- 
suming a form for x^'^^ it) for all so the values of used 
in Eq. (O) were chosen numerically in an ad hoc fashion. 
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FIG. 6. Scaling plot for Pi{l), for the MC simulations 
('MC') and the PCA simulations with Poisson ('PCAP') and 
full ('PCAF') initial conditions. 




FIG. 7. Scaling plot for Pm{m), for the MC simulations 
('MC') and the PCA simulations with Poisson ('PCAP') and 
full ('PCAF') initial conditions. 



III. PROBABILISTIC CELLULAR AUTOMATA 
MODEL 

A. Description of Model 

This model has been described extensively in previous 
publications [^Ul9[| . In the one-dimensional realization of 
this model, there are up to two particles of each species at 
each site, labeled by the direction from which they moved 
onto the site at the previous timestep. The diffusion step 
consists of changing the velocities of these particles, then 
moving the particles onto the neighbouring sites accord- 
ing to their new velocities. If there are two particles per 
site, they both move in opposite directions, whereas a 
single particle will change direction with probability p. 
The value used in these simulations was p = ^, so that 
the particle forgets its previous velocity at each timestep, 
and the model is equivalent to the MC model with Ip = 1. 

The reaction step consists of checking each site for si- 
multaneous occupancy of A and B particles at the start 
of the timestep , and removing any pairs that hopped 
onto the site from opposite directions. Using the segre- 
gated initial condition, and this 'infinite' reaction rate, 
a site can only be occupied by an A-B pair if the A ar- 
rived from the left and the B from the right. This model 
has the same two-sublattice structure as the Monte-Carlo 
model defined above, and this is preserved by the reaction 
algorithm, so these two sublattices must be viewed as two 
independent systems. There is therefore an independent 
nearest-neighbour A-B pair for each sublattice. A multi- 
spin-coding implementation of the algorithm simulates 
64 independent systems at once. 

The quantities Pm{m), Pi{l) and R{x,t) at measure- 
ment time t were estimated by averaging over the interval 
t{l - S) < t < t{l + 6), with S = 0.05. We may estimate 
the order of magnitude of the systematic error that this 



introduces into the measured shape of these quantities. 
Let F{x, t) be the estimate of a function F{x, t) using the 
above method. Then 



1 /•*(l+<5) 

P^T^rJ F{x,t')dt' 



2tS 



t(l-5) 



(19) 



+ lit'-tf^Fix,t) + ..)jdt' (20) 

= F{x,t) + U__F{x,t)+0[{t6f] . (21) 

The fractional error is therefore of order {t5fF/{QF). 

This systematic error has no effect on the scaling be- 
haviour, however. If F{x,t) — t^(t){x/t°'), we have 



F{x,t) 



2tS 



t(l+<5) 



(t')V 



t(l-<5) 



dt' 



(22) 
(23) 



where 0(y) = (2(5)-^ /^^J^/ 6'V(y6'"°)rf«', so F has the 
same scaling properties as F . 

In order to maximize the statistics, the reaction pro- 
file R was measured at every timestep between t(l — i5) 
and t(l -|- 5). However, the quantities m and I are much 
more cumbersome to measure using this program (due to 
the multi-spin coding) , and so were only measured every 
10 timesteps. No significant loss in statistics is incurred, 
since these quantities have very strong time autocorrela- 
tions. The Fortran implementation of this algorithm 
performed 3.7 x 10^ site updates per second on a HP 
9000/715/75 workstation. 
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FIG. 8. Fits of Pm, P;, and R to Eqs. @, (|l|) and ^) 
from the MC simulations ('MC') and PCA simulations with 
Poisson ('PCAF') and full ('PCAF') initial conditions. X-axis 
rescaled and Y-axis shifted for clarity. For curve 'Pmc', see 
text. 

B. Simulation Results 

1. Poisson Initial Condition 

An initial condition with Poisson-like density fluctua- 
tions was prepared by filling each of the appropriate site 
variables (A particle for x < 0, B particles for x > 0) 
with probability j. The lattice size was 4000 sites, and 
at the boundaries particles were free to leave the system, 
with the density at the boundary maintained at an aver- 
age value of i by allowing A particles to enter from the 
left, and B particles to enter from the right, randomly 
with probability j. Measurements were taken at times 
200-102400 timesteps, with the interval between mea- 
surements doubling progressively. The quantities Pm (m) , 
Pm{l) and R{x, t) were measured as described above, and 
then averaged over 82176 independent realizations of the 
system. The quantities m!^''\ l'^'^^ and a;*^''-' were then mea- 
sured from the (l/g)'th power of the normalized q'th mo- 
ment of these quantities. 

Figure ^ shows a plot of {a — b) as a function of {x/t^/'^) 
for three time values, and a plot of T^R-t^/'^ (where S^j = 
/ Rdx) as a function of t. These plots show that, just as 
in the MC simulations, there are no finite size effects. 

Figure |ll| is a log-log plot of ml*-* and I'f'' 
function of time, where 



as a 



(24) 



and fiq is the appropriate scaling factor for Gaussian dis- 
tributions (see Eqs. ( p^ and (1^)). The curves for ml''' 
have been shifted vertically (by 0.2) for clarity, other- 
wise they would be too close to the curves for xl'^' . The 
straight lines are a fit to the last 5 points, for the low- 
est values of q. The gradients of least-square fits for all 
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FIG. 9. Scaling plot of the density difference a{x, t) —b{x, t) 
for the PCA data with random initial condition. Inset: bias 
plot for the total reaction rate E_r = J Rdx. 
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FIG. 10. Scaling plot for R{x,t), for the PCA simulations 
with Poisson ('P') and full ('F') initial conditions. 



the curves are summarized in Table The collapse of 
the curves for different values of q confirms both the 
scalin g hypo thesis and the forms for the scaling func- 
tions (14,1^), and also that the reaction rate profile has 
a Gaussian form: 



R{x,t) 



: exp 



(25) 



Figure ^ shows the effective exponents for x'^' , nrS^^ 
and from the successive gradients in Fig. |l^. The 
curves are much less noisy than those in Fig. ^, by virtue 
of higher statistics and the use of coarse-grained time 
averages. There is a clear trend for the effective exponent 
for x^^' to decrease as time increases, consistent with 
the asymptotic value j predicted elsewhere [|l^,0. The 
exponent for m^^^ appears to increase initially, but the 
last few points appear also to decrease, and in any case 
an asymptotic value 0.375 is ruled out. 

The rescaled forms of P;(Z), Pm(m), and R are denoted 
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FIG. 11. Log-log plot ol ml"', and x'f^ (see text) from 
the PGA simulations with Poisson initial conditions. The 
curves for m'^''-' has been shifted vertically for clarity. 
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FIG. 13. Log-log plot of ^l"', ml"', and a;!"' (see text) from 
the PGA simulations with full initial conditions. The curves 
for m*^"' has been shifted vertically for clarity. 
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FIG. 12. Effective exponents for m'^', and a;'^' (see 
text) from the PGA simulations with Poisson initial condi- 
tions. 
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FIG. 14. Effective exponents for /'^', m'^', and a;^^' (see 
text) from the PGA simulations with full initial conditions. 



by 'PCAP' in Figs. |, @, and respectively. Figure Jsj 
shows a fit of Pi{l), (m) and R to the forms (jl|,p|,|5|^ 
From Eq. (|l9|), using F{x,t) = Af-^ exp{-\jyyt^, the 
fractional error introduced by the coarse-grained time 
averaging is found to be « ((5^/6)(a;/w)'*, where — 
(/ x^Fdx/ J Fdx). The measurement of these quanti- 
ties is therefore expected to be accurate for the first four 
decades or so, as is indeed observed. 



2. Full Initial Conditions 

Because of the exclusion principle in the PCA model, 
the system is completely static in regions where the oc- 
cupation number is zero for one species and assumes its 
maximal value for the other. If one starts from a lattice 
that is filled with A-particles up to x = 0, and filled with 



B-particles for a; > 0, simulations may be speeded up by 
only updating the lattice in the region where a 'hole' has 
penetrated. By checking explicitly that such holes never 
reach the physical boundary of the system, it is possi- 
ble to perform simulations on a system that is effectively 
infinite, so having no finite-size effects. 

Simulations of 64000 independent evolutions of a full 
lattice were run for 409600 timesteps. Measurements of 
Pmi'rn), Pi{l), and R were made using the same method 
as for the Poisson initial condition. Results for these sim- 
ulations are shown in Figs. |l^ and p^. It might be ex- 
pected (considering the arguments in [|16|] ) that this case 
would be in a different universality class from the case 
with randomness in the initial state. However, the results 
for the exponents (see Table |) are very close to those 
measured for the case of Poisson initial conditions, and 
the marked decrease of the exponents for x*^^-' (arguably 
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t5X) 

o 



0.2 




Since the particle profiles vary like x/t^^"^, this means 
that ~ ^1/4, so m ^ i^/s 

In order to assess the validity of this argument, it is 
possible to apply it to a related quantity upon which ana- 
lytical calculations may be made. Consider the diffusion 
equation dtp{x,t) = jd'^p{x,t) in one spatial dimension, 
with an initial condition that consists of a random series 
of negative Dirac delta peaks for a; < and positive Dirac 
delta peaks for a; > 0. That is, 



i=l 



Six + yi), 



(26) 



-3 -2 



-1 



FIG. 15. Plot of a;'^-* and m'^-* against the time-dependent 
current for the two sets of PGA simulations. 



where > 0, j/i > 0. If the intervals between the Xi and 
yi have a Poisson distribution, one has 

(p(a;,0)) =sign(x), (27) 
{pix, 0)p(y, 0)) = sign(xzj) + Six - y), (28) 

where () represents an average over the variables Xi, yi. 
The solution for p may be written in the form 



and are in the 



towards 0.25) is also seen in Fig 
note that the transient trends in m 
opposite sense to the Poisson case. 

Scaling plots for Pi, and R are shown in Figs. ||, 
0, and[Io[ denoted by 'PCAF'. Plots of r^Pi, P„, and 



It is interesting to p(a;,i)= / pix',0)i'Kt) 



\-l/2 



cxp 



ix - x')' 



-1/2 



exp 



(x - Xi)'^ 



R may be found in Fig, 
again have the forms (Q 




confirming that the profiles 
"]!). 

Figure |15| is a plot of m'^^and x^"^^ as a function of 
the time dependent current = J Rdx, from the simu- 
lations both with ('P') and without ('F') Poisson fluc- 
tuations in the initial state. The two curves for a;^^^ 
are almost coincident, which is what would be expected 
if the reaction profile depended upon the current only. 
The curves for m^^\ however, are not quite coincident, 
showing that this quantity is more sensitive to the ini- 
tial condition. Incidentally, numerical tests showed that 
the diffusion current at the origin has Poissonian noise 
whether the initial state contained such fluctuations or 
not. 



IV. THE EFFECT OF POISSONIAN 
FLUCTUATIONS IN THE INITIAL STATE 

The measured value m ^ t^/^ in Ref. was justi- 
fied by an argument about the Poisson fluctuations in 
the initial condition. The argument went as follows: af- 
ter time t, particles within a distance ~ t^/^ have had a 
chance of participating in the reaction. The number of 
particles within a distance t^/'^ is of order t^/'^ ± d}/^. 
Since each reaction event kills precisely one A and one 
B, there is therefore a local surplus ~ t^^^ of one of the 
species. The majority species therefore invades the mi- 
nority species by a distance m, such that the number of 
minority particles between the origin and m is ~ i^/^. 



' exp 



ix + yi) 



2 \ 1 



dx' (29) 



(30) 



Consider the gradient of p 
a,p = ^2(^<)-i/2 [^^exp 



2|p(a:,t) 



-1/2 



t 

x + y, 
t 



Xi exp 



■ exp 



ix - Xi)' 



t 



(31) 



+yi exp - 



ix - x^Y' 



ix + yi 



t 



(32) 



\) is 



The second term on the right hand side of Eq. 
strictly positive. The gradient of p when p is zero is 
therefore strictly positive, so, since p is continuous for all 
t > 0, p is zero at precisely one point, a;o(t) (say). 

It is possible to flnd the probability distribution of 
these zeros, P(a;o), over the ensemble of initial states. 
The position a:o is defined by pixa,t) — 0, or, equiva- 
lently. 



/CO 
p(z, 0) exp 
-oo 



exp 



2za;o 



dz = 0. (33) 



t J \ t 

Suppose that xq ~ i° , where a is expected to be less that 
i, and let e = if'+(i/2)^ ^ith < 6 < (i - a). Then the 
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contribution to the integral in ( |33| ) for \x\ > e is of order 
exp(— e^/i) ~ exp(— i^^), which vanishes as i ^ oo. How- 
ever, for \x\ < e, the argument of the second exponential 
has upper bound 2exo/t — > 0, and so the asymptotic 
value of the integral is found by using the first few terms 
only of the Taylor expansion of this exponential. In other 
words, the leading contribution to xo as i — > oo is given 

by 



p{z,0)e--dz / p{z,0)ze- — dz=^0. (34) 

) ^ J — oo 

The expectation value of the second moment of Xq is 



1.25 



[Xn) = — 



/ J p{x, 0)e dx J p{y, 0)e "* 



Jxpix,0)e dx J yp{y,0)e ^dy, 



(35) 



To evaluate this average, write p[x, 0) — sign(a::) + t{x), 
where {t{x)) = and {T{x)T{y)) = S{x — y). Then 
/ xp{x,Q)exp{—x'^/t)dx = t + J xt{x) exp{—x^ /t)dx, the 
second term being typically much smaller than the first. 
To find the leading contribution to xq, it is sufficient to 
replace J xp{x, 0) exp(— a;^/i)(ix in the denominator by t. 
We therefore have 



oc />oo 



(p(x,0)p(y,0)) 



X exp 



' dx 



x^ + y^ 



dxdy 



1 irt 

Similarly, the 27z'th moment of P(xo) is of the form 

(^o"> ^ ^ dxi... dX2n 

2" 



(36) 
(37) 

(38) 
(39) 



i=l 



X exp 

\ I 

X (p(xi,0)...p(x2„,0)) + ... (40) 

22»n!^8; ^ ' 

For a dis- 

tribution of the form P(a;o) = \/ (A/tt) exp(— Ax§), one 
has 



„2n\ 



(2n)! 

Comparison with ([4l|) gives 



(42) 




FIG. 16. Log-log plot of {a;^) versus t for q = 2 { 
8 (A) and 16 (x), from numerical solutions of Eq. (I 
straight lines are the asymptotic solutions from Eq, 



^(-o) = ^/fexp [-{^^xl^. 




(43) 



The distribution of xq is therefore characterized by a sin- 
gle lengthscale A""'^/^ oc i^/^. 

Figure |l^ shows the moments of Poi averaged over 
10000 realizations, from a numerical solution of the zero 
of p from Eq. (pO), compared with the asymptotic pre- 
dictions of Eq. ( |4lD . 

From Eq. (|l|), the density difference (a — b) in the 
reaction-diffusion problem, averaged over evolutions, 
obeys a simple diffusion equation. The quantity p with 
the initial condition (|2^) is therefore equal to (a — b) for 
the initial condition with Poisson fluctuations used in the 
numerical simulations, with negative peaks correspond- 
ing to A particles and positive peaks corresponding to B 
particles. The quantity differs from ■m{t) because the 
latter contains further fluctuations due to the diffusive 
noise that has been averaged over in the former. How- 
ever, the argument used in ||l^ to obtain m ~ t^^^ may 
be applied equally well to xq. The reaction centre shifts 
to compensate for a local majority of order t^/^ in one of 
the species, and the argument predicts xq ~ t^^^. It is 
interesting to note that the correct exponent is obtained 
if the initial value a{x, 0) — oq is used instead of the 
value a{x^t) oc x/t^^'^ at time t in the balance equation 
J^" a{x, t) ~ t^/^. This ambiguity is probably the reason 
for the argument being incorrect. 



V. CONCLUSIONS 

It appears from extensive simulations that the reac- 
tion profile in this system has the same simple dynamic 
scaling form, independently of the presence of an exclu- 
sion principle and of randomness in the initial state. The 
motion of the reaction centre due to the Poisson noise ap- 
pears only to account for a contribution of order ^ t^/* 
to the reaction width, which is not large enough to al- 
ter the scaling behaviour. The measured exponent ~ 0.29 
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TABLE I. Comparison of the simulation results in this pa- 
per for the Monte-Carlo model (MC) and PCA model with 
Poisson (PCAP) and Full (PCAF) initial conditions with 
those of Araujo et al (ALHS) [^. Numbers in parentheses 
represent the statistical error in the preceding digit. 





ALHS 


MC 


PCAP 


PCAF 


Size 


2000 


4000 


4000 


Infinite 


Exclusion prle.? 


No 


No'' 


Yes 


Yes 


Initial density 


1.0 


1.0 


0.5 


2.0 


Initial State 


(Uniform) 


Poisson 


Poisson 


Uniform 


Averaging 


dOOO-IoOUU 


21000 


82176 


64000 


Max time 


20000 


zoOOO 


1 AO A AA 

102400 


409600 


Exponents: 












0.25 


0.251(3) 


0.2510(6) 


0.2542(4) 


;(16) 


0.25 


0.23(1) 


0.248(3) 


0.2609(3) 




0.375 


0.281(4) 


0.287(1) 


0.300(1) 




0.375 


0.29(1) 


0.284(1) 


0.299(3) 




0.312 


0.2799(2)'" 


0.286(2) 


0.291(1) 




0.359 


0.282(2)*" 


0.280(4) 


0.293(1) 




0.367 


0.30(2)"" 


0.28(1) 


0.293(2) 


Fit over last. . . 


N/A 


8 points 


5 points 


6 points 



they used would therefore not appear to be vahd, unless 
it describes a regime inaccessible to simulations. 

The simulation evidence in favour of dynamic scaling 
in this model is very strong. However, the numerical 
evidence that all lengthscales scale asymptotically as t^/* 
is far from conclusive, and so needs to be put on a sound 
theoretical basis, either by an exact calculation or by a 
rigorous justification for the analogy with the static case. 
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''See text 

""Measured from X*"' 



describing both the reaction width and the midpoint fluc- 
tuations appears in fact to be decreasing slowly in time, 
with favourable evidence for an asymptotic value 0.25. 
This, together with the measured form for the reaction 
profile, is consistent with the steady-state results being 
applicable |l|,0j2O|. 

It is, nevertheless, surprising that the approach to the 
asymptotic behaviour should be so slow. It is not clear 
whether logarithmic corrections should be present, as 
they do not occur in the steady-state problem . How- 
ever, in these simulations the ratio of the reaction width 
w to the diffusion length (£)/:) was never smaller than 
K, 0.2, whereas the application of the steady-state argu- 
ment requires that this ratio be small. This could ac- 
count for the fact that the asymptotic regime has not 
been reached. Simulations where this ratio is truly small 
would not appear to be practical at present. 

An investigation of the simulation procedure used by 
Araujo et al has revealed a few errors in the results pub- 
lished in p6[ | . A repeat of their simulations appears to 
confirm the results of the present article for Pm and P/, 
and the behaviour m ~ t'^-^°, but does not find that R 
satisfies a scaling ansatz |2l|| . This inconsistency between 
my results and those of Araujo et al is currently unex- 
plained. 



A recent calculation by Rodriguez and Wio |22| sug- 
gests that the reaction profile R should be the superpo- 
sition of two scaling forms, with width exponents ^ and 
I respectively. However, these exponents and the form 
they predict for R (~ exp[— (x/w)'^/^]) do not agree with 
the results of simulations. The approximation scheme 
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